week 8: multilevel models
MELSM
double-advanced potions multilevel models
We’ll use data provided by Williams and colleagues (2021 .
<- read_csv ("https://raw.githubusercontent.com/sjweston/uobayes/refs/heads/main/files/data/external_data/williams.csv" )# scaled time variable <- d %>% mutate (day01 = (day - 2 ) / max ((day - 2 )))distinct (d, id) %>% count ()
# A tibble: 1 × 1
n
<int>
1 193
%>% count (id) %>% summarise (median = median (n),min = min (n),max = max (n))
# A tibble: 1 × 3
median min max
<int> <int> <int>
1 74 8 99
Code
%>% count (id) %>% ggplot (aes (x = n)) + geom_histogram (fill = "#1c5253" , color = "white" ) + scale_x_continuous ("number of days" , limits = c (0 , NA ))
mean sd 5.5% 94.5% histogram
id 98.61029694 63.7493291 10.00000000 207.0000000 ▇▇▇▃▅▅▅▃▃▃▂▂
female 0.57016803 0.4950710 0.00000000 1.0000000 ▅▁▁▁▁▁▁▁▁▇
PA.std 0.01438236 1.0241384 -1.68129708 1.7514660 ▁▁▃▇▇▃▁▁
day 44.17962096 27.6985612 6.00000000 92.0000000 ▇▇▇▇▅▅▅▃▃▃
PA_lag 0.01992587 1.0183833 -1.72043510 1.7170361 ▁▂▃▅▇▇▅▃▂▁
NA_lag -0.04575229 0.9833161 -0.87507181 1.9904675 ▇▃▂▁▁▁▁▁▁▁▁▁▁▁
steps.pm 0.05424387 0.6298941 -1.02580678 1.0113559 ▁▂▇▇▃▁▁▁
steps.pmd 0.02839160 0.7575395 -1.12359512 1.2299739 ▁▁▇▇▁▁▁▁▁▁▁▁▁▁
NA.std -0.07545093 0.9495660 -1.04842926 1.8110607 ▁▁▁▇▂▁▁▁▁▁
day01 0.43040430 0.2826384 0.04081633 0.9183673 ▇▇▇▇▅▅▅▃▃▃
Code
set.seed (14 )%>% nest (data= - id) %>% slice_sample (n = 6 ) %>% unnest (data) %>% ggplot (aes (x = day, y = NA_lag)) + geom_line (color = "grey" ) + geom_point (color = "#1c5253" , size = 1 / 2 ) + geom_line (aes (y = PA_lag), color = "darkgrey" ) + geom_point (aes (y = PA_lag), color = "#e07a5f" , size = 1 / 2 ) + ylab ("affect (standardized)" ) + facet_wrap (~ id)
Let’s fit an MLM with varying intercepts and slopes.
Likelihood function and linear model
\[\begin{align*}
\text{NA}_i &\sim \text{Normal}(\mu_i, \sigma) \\
\mu_i &= \alpha_{\text{id}[i]} + \beta_{1\text{id}[i]}\text{NA_lag}_i + \beta_{2\text{id}[i]}\text{PA_lag}_i
\end{align*}\]
Varying intercepts and slopes:
\[\begin{align*}
\alpha_{\text{id}[i]} &= \alpha + u_{\alpha,\text{id}[i]} \\
\beta_{1\text{id}[i]} &= \beta + u_{\beta,\text{id}[i]} \\
\beta_{2\text{id}[i]} &= \beta + u_{\beta,\text{id}[i]}
\end{align*}\]
Random effects:
\[\begin{align*}
\begin{bmatrix} u_{\alpha,\text{id}[i]} \\ u_{\beta,1\text{id}[i]} \\ u_{\beta,2\text{id}[i]} \end{bmatrix} &\sim \text{MVNormal}\left(\begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix}, \mathbf{S}\right) \\
\mathbf{S} &= \begin{pmatrix} \sigma_{\alpha} & 0 \\ 0 \\ 0 & \sigma_{\beta}\end{pmatrix}\mathbf{R}\begin{pmatrix} \sigma_{\alpha} & 0 \\ 0 \\ 0 & \sigma_{\beta}\end{pmatrix}
\end{align*}\]
Priors: \[\begin{align*}
\alpha &\sim \text{Normal}(0,0.2) \\
\beta &\sim \text{Normal}(0,0.5) \\
\sigma &\sim \text{Exponential}(1) \\
\sigma_{\alpha} &\sim \text{Exponential}(1) \\
\sigma_{\beta} &\sim \text{Exponential}(1) \\
\mathbf{R} &\sim \text{LKJcorr}(2)
\end{align*}\]
<- brm (data = d,family = gaussian,~ 1 + NA_lag + PA_lag + (1 + NA_lag + PA_lag | id),prior = c (prior (normal (0 , 0.2 ), class = Intercept),prior (normal (0 , 1 ), class = b),prior (exponential (1 ), class = sd),prior (exponential (1 ), class = sigma),prior (lkj (2 ), class = cor)),iter = 5000 , warmup = 1000 , chains = 4 , cores = 4 ,seed = 14 ,file = here ("files/models/m82.1" ))
What if I wanted to model positive affect? Usually, I’d have to fit a second model.
Multivariate MLMs
But what if I told you I could estimate both positive and negative affect simultaneously?
Likelihood function
\[\begin{align*}
\text{NA}_i &\sim \text{Normal}(\mu_{\text{NA},i}, \sigma_{\text{NA}}) \\
\mu_{\text{NA},i} &= \alpha_{\text{NA},\text{id}[i]} + \beta_{1\text{NA},\text{id}[i]}\,\text{PA_lag}_i + \beta_{2\text{NA},\text{id}[i]}\,\text{NA_lag}_i \\
\\
\text{PA}_i &\sim \text{Normal}(\mu_{\text{PA},i}, \sigma_{\text{PA}}) \\
\mu_{\text{PA},i} &= \alpha_{\text{PA},\text{id}[i]} + \beta_{1\text{PA},\text{id}[i]}\,\text{PA_lag}_i + \beta_{2\text{PA},\text{id}[i]}\,\text{NA_lag}_i
\end{align*}\]
Varying intercepts and slopes:
\[\begin{align*}
\alpha_{\text{NA},\text{id}[i]} &= \alpha_{\text{NA}} + u_{\alpha,1\text{NA},\text{id}[i]} \\
\beta_{1\text{NA},\text{id}[i]} &= \beta_{1\text{NA}} + u_{\beta,1\text{NA},\text{id}[i]} \\
\beta_{2\text{NA},\text{id}[i]} &= \beta_{2\text{NA}} + u_{\beta,2\text{NA},\text{id}[i]} \\
\alpha_{\text{PA},\text{id}[i]} &= \alpha_{\text{PA}} + u_{\alpha,\text{PA},\text{id}[i]} \\
\beta_{1\text{PA},\text{id}[i]} &= \beta_{1\text{PA}} + u_{1\beta,\text{PA},\text{id}[i]} \\
\beta_{2\text{PA},\text{id}[i]} &= \beta_{2\text{PA}} + u_{2\beta,\text{PA},\text{id}[i]}
\end{align*}\]
Random Effects: Multivariate Distribution
\[\begin{align*}
\mathbf{u}_{\text{id}[i]} =
\begin{bmatrix}
u_{\alpha,\text{NA},\text{id}[i]} \\
u_{\beta,1\text{NA},\text{id}[i]} \\
u_{\beta,2\text{NA},\text{id}[i]} \\
u_{\alpha,\text{PA},\text{id}[i]} \\
u_{\beta,1\text{PA},\text{id}[i]} \\
u_{\beta,2\text{PA},\text{id}[i]}
\end{bmatrix}
&\sim \text{MVNormal}(\mathbf{0}, \boldsymbol{\Sigma}) \\
\boldsymbol{\Sigma} &= \mathbf{L} \mathbf{R} \mathbf{L} \quad
\\ \text{with} \quad
\mathbf{L} &= \text{diag}(\sigma_{\alpha,\text{NA}}, \sigma_{\beta,1\text{NA}}, \sigma_{\beta,2\text{NA}}, \sigma_{\alpha,\text{PA}}, \sigma_{\beta,1\text{PA}}, \sigma_{\beta,2\text{PA}})
\end{align*}\]
Priors
\[\begin{align*}
\alpha_{\text{NA}}, \alpha_{\text{PA}} &\sim \text{Normal}(0, 0.2) \\
\beta_{1\text{NA}}, \beta_{2\text{NA}}, \beta_{1\text{PA}}, \beta_{2\text{PA}} &\sim \text{Normal}(0, 0.5) \\
\sigma_{\text{NA}}, \sigma_{\text{PA}} &\sim \text{Exponential}(1) \\
\sigma_{\alpha,\text{NA}}, \sigma_{\beta,1\text{NA}}, \sigma_{\beta,2\text{NA}}, \sigma_{\alpha,\text{PA}}, \sigma_{\beta,1\text{PA}}, \sigma_{\beta,2\text{PA}} &\sim \text{Exponential}(1) \\
\mathbf{R} &\sim \text{LKJcorr}(2)
\end{align*}\]
= bf (NA.std ~ 1 + NA_lag + PA_lag + (1 + NA_lag + PA_lag | c | id))= bf (PA.std ~ 1 + NA_lag + PA_lag + (1 + NA_lag + PA_lag | c | id))<- brm (data = d,family = gaussian,+ bf_pa + set_rescor (TRUE ),prior = c (prior (normal (0 , 0.2 ), class = Intercept),prior (normal (0 , 1 ), class = b),prior (lkj (2 ), class = cor)),iter = 5000 , warmup = 1000 , chains = 4 , cores = 4 ,seed = 14 ,file = here ("files/models/m82.2" ))
Family: MV(gaussian, gaussian)
Links: mu = identity; sigma = identity
mu = identity; sigma = identity
Formula: NA.std ~ 1 + NA_lag + PA_lag + (1 + NA_lag + PA_lag | c | id)
PA.std ~ 1 + NA_lag + PA_lag + (1 + NA_lag + PA_lag | c | id)
Data: d (Number of observations: 13033)
Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 1;
total post-warmup draws = 16000
Multilevel Hyperparameters:
~id (Number of levels: 193)
Estimate Est.Error l-95% CI u-95% CI Rhat
sd(NAstd_Intercept) 0.44 0.03 0.39 0.50 1.00
sd(NAstd_NA_lag) 0.15 0.01 0.12 0.18 1.00
sd(NAstd_PA_lag) 0.07 0.02 0.03 0.10 1.00
sd(PAstd_Intercept) 0.43 0.03 0.38 0.48 1.00
sd(PAstd_NA_lag) 0.05 0.02 0.01 0.08 1.00
sd(PAstd_PA_lag) 0.11 0.01 0.09 0.14 1.00
cor(NAstd_Intercept,NAstd_NA_lag) 0.09 0.10 -0.10 0.29 1.00
cor(NAstd_Intercept,NAstd_PA_lag) 0.30 0.15 -0.00 0.60 1.00
cor(NAstd_NA_lag,NAstd_PA_lag) 0.07 0.18 -0.30 0.40 1.00
cor(NAstd_Intercept,PAstd_Intercept) -0.13 0.08 -0.28 0.02 1.00
cor(NAstd_NA_lag,PAstd_Intercept) -0.00 0.11 -0.22 0.22 1.00
cor(NAstd_PA_lag,PAstd_Intercept) -0.12 0.16 -0.42 0.18 1.01
cor(NAstd_Intercept,PAstd_NA_lag) -0.06 0.17 -0.39 0.28 1.00
cor(NAstd_NA_lag,PAstd_NA_lag) -0.03 0.21 -0.43 0.38 1.00
cor(NAstd_PA_lag,PAstd_NA_lag) 0.26 0.24 -0.24 0.69 1.00
cor(PAstd_Intercept,PAstd_NA_lag) -0.24 0.18 -0.59 0.12 1.00
cor(NAstd_Intercept,PAstd_PA_lag) -0.18 0.11 -0.39 0.04 1.00
cor(NAstd_NA_lag,PAstd_PA_lag) 0.10 0.15 -0.19 0.39 1.00
cor(NAstd_PA_lag,PAstd_PA_lag) -0.12 0.18 -0.48 0.24 1.00
cor(PAstd_Intercept,PAstd_PA_lag) 0.09 0.11 -0.11 0.30 1.00
cor(PAstd_NA_lag,PAstd_PA_lag) 0.32 0.21 -0.13 0.68 1.00
Bulk_ESS Tail_ESS
sd(NAstd_Intercept) 3584 6292
sd(NAstd_NA_lag) 6665 10373
sd(NAstd_PA_lag) 1557 1372
sd(PAstd_Intercept) 4387 7836
sd(PAstd_NA_lag) 2112 2250
sd(PAstd_PA_lag) 7127 10096
cor(NAstd_Intercept,NAstd_NA_lag) 6685 9729
cor(NAstd_Intercept,NAstd_PA_lag) 7006 7448
cor(NAstd_NA_lag,NAstd_PA_lag) 7576 7022
cor(NAstd_Intercept,PAstd_Intercept) 2349 4989
cor(NAstd_NA_lag,PAstd_Intercept) 897 1718
cor(NAstd_PA_lag,PAstd_Intercept) 280 469
cor(NAstd_Intercept,PAstd_NA_lag) 11866 10018
cor(NAstd_NA_lag,PAstd_NA_lag) 8743 9256
cor(NAstd_PA_lag,PAstd_NA_lag) 2757 5270
cor(PAstd_Intercept,PAstd_NA_lag) 14554 10682
cor(NAstd_Intercept,PAstd_PA_lag) 12167 12720
cor(NAstd_NA_lag,PAstd_PA_lag) 4607 8122
cor(NAstd_PA_lag,PAstd_PA_lag) 2163 3286
cor(PAstd_Intercept,PAstd_PA_lag) 14121 12777
cor(PAstd_NA_lag,PAstd_PA_lag) 1588 2031
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
NAstd_Intercept -0.03 0.03 -0.10 0.03 1.00 1348 2892
PAstd_Intercept -0.01 0.03 -0.07 0.05 1.00 2691 5026
NAstd_NA_lag 0.35 0.02 0.32 0.38 1.00 7233 10519
NAstd_PA_lag 0.01 0.01 -0.01 0.03 1.00 7186 11061
PAstd_NA_lag 0.04 0.01 0.02 0.06 1.00 16764 13101
PAstd_PA_lag 0.45 0.01 0.43 0.48 1.00 12268 12423
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma_NAstd 0.59 0.00 0.59 0.60 1.00 27696 10725
sigma_PAstd 0.54 0.00 0.54 0.55 1.00 31518 11100
Residual Correlations:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
rescor(NAstd,PAstd) -0.02 0.01 -0.04 -0.01 1.00 32445 12034
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
These models can help us address many questions. First, let’s decompose the between- and within-person variance of these outcomes.
Code
# Compute person-level means and residuals <- d %>% group_by (id) %>% mutate (N_A_mean = mean (NA.std, na.rm = TRUE ),P_A_mean = mean (PA.std, na.rm = TRUE ),N_A_within = NA.std - N_A_mean,P_A_within = PA.std - P_A_mean%>% ungroup ()<- d_resid %>% :: select (id, N_A_mean, P_A_mean) %>% distinct ()<- cor (d$ NA.std, d$ PA.std, use= "complete.obs" )<- cor (between_df$ N_A_mean, between_df$ P_A_mean, use = "complete.obs" )<- cor (d_resid$ N_A_within, d_resid$ P_A_within, use = "complete.obs" )<- m2 %>% spread_draws (rescor__NAstd__PAstd) %>% mean_qitibble (Level = c ("Total" , "Between-person" , "Within-person" , "Residual" ),Correlation = c (total_cor, between_cor, within_cor, as.numeric (residual_cor[1 ,1 ])))
# A tibble: 4 × 2
Level Correlation
<chr> <dbl>
1 Total -0.0597
2 Between-person -0.0544
3 Within-person -0.0460
4 Residual -0.0249
Code
%>% gather_draws (r_id__NAstd[id, term], r_id__PAstd[id, term]) %>% mutate (.variable= str_remove (.variable, "r_id__" )) %>% pivot_wider (names_from = .variable, values_from = .value) %>% # just for presenting to class arrange (.draw, id, term) %>% select (- .chain, - .iteration) %>% mutate (across (c (NAstd, PAstd), \(x) round (x, 2 )))
# A tibble: 9,264,000 × 5
# Groups: id, term [579]
id term .draw NAstd PAstd
<int> <chr> <int> <dbl> <dbl>
1 1 Intercept 1 0.14 0.05
2 1 NA_lag 1 0.11 0.08
3 1 PA_lag 1 0.07 0.13
4 2 Intercept 1 -0.4 0.88
5 2 NA_lag 1 -0.03 0
6 2 PA_lag 1 -0.04 0.16
7 3 Intercept 1 0.82 0
8 3 NA_lag 1 0.33 -0.09
9 3 PA_lag 1 0.1 -0.2
10 4 Intercept 1 0.2 -0.08
# ℹ 9,263,990 more rows
spaghetti plots
Code
= data.frame (PA_lag = seq (from = min (d$ PA_lag), to= max (d$ PA_lag), length.out= 50 ),NA_lag = 0 ,id = max (d$ id) + 1 %>% add_epred_draws (m2, allow_new_levels= T) %>% filter (.draw <= 50 ) %>% ggplot (aes (x = PA_lag, y = .epred)) + geom_line (aes (group = .draw, color = .category), alpha= .2 ) + scale_color_manual (values = c ("#1c5253" , "#e07a5f" )) + facet_wrap (~ .category) + guides (color= "none" ) + labs (x = "Positive Affect (lagged)" , y = "Expected score" )
Code
= data.frame (NA_lag = seq (from = min (d$ PA_lag), to= max (d$ PA_lag), length.out= 50 ),PA_lag = 0 ,id = max (d$ id) + 1 %>% add_epred_draws (m2, allow_new_levels= T) %>% filter (.draw <= 50 ) %>% ggplot (aes (x = NA_lag, y = .epred)) + geom_line (aes (group = .draw, color = .category), alpha= .2 ) + scale_color_manual (values = c ("#1c5253" , "#e07a5f" )) + facet_wrap (~ .category) + guides (color= "none" ) + labs (x = "Negative Affect (lagged)" , y = "Expected score" )
But of course, we’ve allowed these effects to vary for each individual.
= expand.grid (PA_lag = seq (from = min (d$ PA_lag), to= max (d$ PA_lag), length.out= 50 ),NA_lag = 0 ,id = sample (d$ id, size = 4 , replace= F)%>% add_epred_draws (m2) %>% filter (.draw <= 20 ) %>% ggplot (aes (x = PA_lag, y = .epred)) + geom_line (aes (group = .draw, color = .category), alpha= .2 ) + scale_color_manual (values = c ("#1c5253" , "#e07a5f" )) + facet_grid (.category~ id) + guides (color= "none" ) + labs (x = "Positive Affect (lagged)" , y = "Expected score" )
We may even want to see whether there’s a correlation between the effect of positive affect on PA and on NA.
= m2 %>% spread_draws (r_id__PAstd[id, term],%>% filter (term == "PA_lag" ) %>% %>% ggplot (aes (x = r_id__PAstd, y = r_id__NAstd)) + geom_point () + labs (x = "Effect of PA on PA" , y= "Effect of PA on NA" )
But of course, each of these points is a summary of a distribution! The real benefit of modeling these outcomes jointly is that we can also get a distribution of the correlations between these effects.
%>% spread_draws (` cor.* ` , regex = T) %>% pivot_longer (starts_with ("cor" ),names_prefix = "cor_id__" ) %>% group_by (name) %>% mean_qi (value) %>% arrange (desc (abs (value)))
# A tibble: 15 × 7
name value .lower .upper .width .point .interval
<chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 PAstd_NA_lag__PAstd_PA_lag 0.318 -0.132 0.681 0.95 mean qi
2 NAstd_Intercept__NAstd_PA_l… 0.303 -0.00369 0.601 0.95 mean qi
3 NAstd_PA_lag__PAstd_NA_lag 0.261 -0.240 0.692 0.95 mean qi
4 PAstd_Intercept__PAstd_NA_l… -0.241 -0.588 0.123 0.95 mean qi
5 NAstd_Intercept__PAstd_PA_l… -0.183 -0.391 0.0356 0.95 mean qi
6 NAstd_Intercept__PAstd_Inte… -0.133 -0.284 0.0217 0.95 mean qi
7 NAstd_PA_lag__PAstd_PA_lag -0.123 -0.479 0.242 0.95 mean qi
8 NAstd_PA_lag__PAstd_Interce… -0.121 -0.419 0.180 0.95 mean qi
9 NAstd_NA_lag__PAstd_PA_lag 0.102 -0.190 0.393 0.95 mean qi
10 PAstd_Intercept__PAstd_PA_l… 0.0950 -0.113 0.296 0.95 mean qi
11 NAstd_Intercept__NAstd_NA_l… 0.0946 -0.103 0.291 0.95 mean qi
12 NAstd_NA_lag__NAstd_PA_lag 0.0685 -0.300 0.405 0.95 mean qi
13 NAstd_Intercept__PAstd_NA_l… -0.0645 -0.393 0.283 0.95 mean qi
14 NAstd_NA_lag__PAstd_NA_lag -0.0334 -0.433 0.383 0.95 mean qi
15 NAstd_NA_lag__PAstd_Interce… -0.00442 -0.223 0.220 0.95 mean qi
modeling variability
Traditionally, psycho-social research has focused on identifying trait-like behaviors - for example, how negative a person is on average. But more recent work has broadened to include states - that is, whether and how much individuals fluctuate in their thoughts, feelings, and behaviors over time.
On the macro time scale, we typically assume stable traits, but on the micro time scale (day to day), we often observe substantial variability in these same traits. This presents a methodological challenge: how do we simultaneously model both the stable personality traits (the “location”) and the fluctuations (the “scale”)? The standard approach has been a two-stage process:
compute individual means (iMs) in a mixed effects model
obtain individual standard deviations (iSDs) from the residuals
But this approach has several drawbacks:
It can result in unreliable estimates
It’s particularly sensitive to the number of measurement occasions
It assumes independence of means and variances, which seems unlikely
It results in biased variance estimates
mixed-effects location scale models
In a standard linear mixed effects model with repeated measurements, we have:
\[
y_i = X_i\beta + Z_ib_i + \epsilon_i
\] Where:
\(y_i\) is the observations for the ith person\(X_i\beta\) represents the fixed effects\(Z_ib_i\) represents the random effects (person-specific deviations)\(\epsilon_i\) represents the error term (or states)
The key limitation of this standard model is that it assumes the same error variance \((\sigma^2)\) for everyone. The MELSM, instead, allows \(\sigma^2\) to differ at the individual level and even at different time points. The scale component is modeled as:
\[
\varphi_i = \text{exp}(W_i\eta + V_it_i)
\] Where:
\(\varphi_i\) contains the error variances\(\eta\) represents fixed effects for the scale\(t_i\) represents random effects for the scaleThe exponential function ensures positive variance estimates
Let’s go back to a model with one outcome, but we can make it even better.
\[\begin{align*}
\text{NA}_{ij} &\sim \text{Normal}(\mu_{ij}, \sigma) \\
\mu_{ij} &= \beta_0 + \beta_1\text{time}_{ij} + \mu_{0i} + \mu_{1i}\text{time}_{ij}\\
\text{log}(\sigma_i) &= \eta_0 + \mu_{2i}
...
\end{align*}\]
<- brm (data = d,family = gaussian,bf (NA.std ~ 1 + day01 + (1 + day01 | i| id),~ 1 + (1 | i| id)),prior = c (prior (normal (0 , 0.2 ), class = Intercept),prior (normal (0 , 1 ), class = b),prior (exponential (1 ), class = sd),prior (normal (0 ,1 ), class = Intercept, dpar= sigma),prior (exponential (1 ), class = sd, dpar= sigma),prior (lkj (2 ), class = cor)),iter = 5000 , warmup = 1000 , chains = 4 , cores = 4 ,seed = 14 ,file = here ("files/models/m82.3" ))
We should note a few things about the brm() syntax. First, because we modeled both \(\mu_{ij}\) and \(\sigma_i\) , we nested both model formulas within the bf() function. Second, because the brms default is to use the log link when modeling \(\sigma_i\) , there was no need to explicitly set it that way in the family line. However, we could have if we wanted to. Third, notice our use of the |i| syntax within the parentheses in the formula lines. If we had used the conventional | syntax, that would have not allowed our \(\mu_{ij}\) parameters to correlate with\(\mu_{0i}\) and \(\mu_{1i}\) from the mean structure. It would have effectively set . Finally, notice how within the prior() functions, we explicitly referred to those for the new \(\sigma\) structure with the dpar = sigma operator.
